{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import os\n",
    "import math\n",
    "import time\n",
    "import yaml\n",
    "import pynvml\n",
    "os.environ['CUDA_VISIBLE_DEVICES'] = '5'\n",
    "import sys\n",
    "import random\n",
    "sys.path.append('/home/yang_liu/python_workspace/3DGS')\n",
    "import torchvision\n",
    "from os import makedirs\n",
    "from torch import nn\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm\n",
    "from scene.datasets import GSDataset\n",
    "from scene import LargeScene, GaussianModel, GaussianModelLOD\n",
    "from gaussian_renderer import render, render_v2, render_v3\n",
    "from diff_gaussian_rasterization import GaussianRasterizationSettings, GaussianRasterizer\n",
    "from arguments import GroupParams\n",
    "from utils.general_utils import PILtoTorch\n",
    "from utils.graphics_utils import getWorld2View2, getProjectionMatrix\n",
    "from utils.large_utils import which_block, block_filtering\n",
    "from utils.sh_utils import eval_sh\n",
    "from transforms3d.euler import euler2mat, mat2euler\n",
    "\n",
    "WARNED = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BlockedGaussian:\n",
    "\n",
    "    gaussians : GaussianModel\n",
    "\n",
    "    def __init__(self, gaussians, lp, range=[0, 1], scale=1.0):\n",
    "        self.xyz = []\n",
    "        self.opacity = []\n",
    "        self.rotation = []\n",
    "        self.scaling = []\n",
    "        self.features = []\n",
    "        self.cov3D = []\n",
    "        self.cell_corners = []\n",
    "        self.max_sh_degree = lp.sh_degree\n",
    "        self.device = gaussians.get_xyz.device\n",
    "\n",
    "        self.block_dim = lp.block_dim\n",
    "        self.num_cell = lp.block_dim[0] * lp.block_dim[1] * lp.block_dim[2]\n",
    "        self.aabb = lp.aabb\n",
    "        self.scale = scale\n",
    "        self.range = range\n",
    "\n",
    "        self.cell_divider(gaussians)\n",
    "        self.cell_corners = torch.stack(self.cell_corners, dim=0)\n",
    "\n",
    "    def cell_divider(self, gaussians):\n",
    "        with torch.no_grad():\n",
    "            xyz = gaussians.get_xyz\n",
    "            for cell_idx in range(self.num_cell):\n",
    "                cell_mask = block_filtering(cell_idx, xyz, self.aabb, self.block_dim, self.scale)\n",
    "                self.xyz.append(xyz[cell_mask])\n",
    "                self.opacity.append(gaussians.get_opacity[cell_mask])\n",
    "                self.rotation.append(gaussians.get_rotation[cell_mask])\n",
    "                self.scaling.append(gaussians.get_scaling[cell_mask])\n",
    "                self.features.append(gaussians.get_features[cell_mask])\n",
    "                self.cov3D.append(gaussians.get_covariance(self.scale)[cell_mask])\n",
    "                xyz_min = xyz[cell_mask].min(dim=0)[0]\n",
    "                xyz_max = xyz[cell_mask].max(dim=0)[0]\n",
    "                corners = torch.tensor([[xyz_min[0], xyz_min[1], xyz_min[2]],\n",
    "                                       [xyz_min[0], xyz_min[1], xyz_max[2]],\n",
    "                                       [xyz_min[0], xyz_max[1], xyz_min[2]],\n",
    "                                       [xyz_min[0], xyz_max[1], xyz_max[2]],\n",
    "                                       [xyz_max[0], xyz_min[1], xyz_min[2]],\n",
    "                                       [xyz_max[0], xyz_min[1], xyz_max[2]],\n",
    "                                       [xyz_max[0], xyz_max[1], xyz_min[2]],\n",
    "                                       [xyz_max[0], xyz_max[1], xyz_max[2]]], device=xyz.device)\n",
    "                self.cell_corners.append(corners)\n",
    "    \n",
    "    def get_xyz(self, indices):\n",
    "        out = torch.tensor([], device=self.device)\n",
    "        if len(indices) > self.range[0]:\n",
    "            indice_end = min(len(indices), self.range[1])\n",
    "            out = torch.cat([self.xyz[indices[i].item()] for i in range(self.range[0], indice_end)])\n",
    "        return out\n",
    "\n",
    "    def get_opacity(self, indices):\n",
    "        out = torch.tensor([], device=self.device)\n",
    "        if len(indices) > self.range[0]:             \n",
    "            indice_end = min(len(indices), self.range[1])\n",
    "            out = torch.cat([self.opacity[indices[i].item()] for i in range(self.range[0], indice_end)])\n",
    "        return out\n",
    "    \n",
    "    def get_rotation(self, indices):\n",
    "        out = torch.tensor([], device=self.device)\n",
    "        if len(indices) > self.range[0]:             \n",
    "            indice_end = min(len(indices), self.range[1])\n",
    "            out = torch.cat([self.rotation[indices[i].item()] for i in range(self.range[0], indice_end)])\n",
    "        return out\n",
    "    \n",
    "    def get_scaling(self, indices):\n",
    "        out = torch.tensor([], device=self.device)\n",
    "        if len(indices) > self.range[0]:             \n",
    "            indice_end = min(len(indices), self.range[1])\n",
    "            out = torch.cat([self.scaling[indices[i].item()] for i in range(self.range[0], indice_end)])\n",
    "        return out\n",
    "    \n",
    "    def get_features(self, indices):\n",
    "        out = torch.tensor([], device=self.device)\n",
    "        if len(indices) > self.range[0]:             \n",
    "            indice_end = min(len(indices), self.range[1])\n",
    "            out = torch.cat([self.features[indices[i].item()] for i in range(self.range[0], indice_end)])\n",
    "        return out\n",
    "    \n",
    "    def get_cov3D(self, indices):\n",
    "        out = torch.tensor([], device=self.device)\n",
    "        if len(indices) > self.range[0]:             \n",
    "            indice_end = min(len(indices), self.range[1])\n",
    "            out = torch.cat([self.cov3D[indices[i].item()] for i in range(self.range[0], indice_end)])\n",
    "        return out\n",
    "\n",
    "class BlockedGaussianV2:\n",
    "\n",
    "    gaussians : GaussianModel\n",
    "\n",
    "    def __init__(self, gaussians, lp, range=[0, 1], scale=1.0, compute_cov3D_python=False):\n",
    "        self.feats = []\n",
    "        self.cell_corners = []\n",
    "        self.max_sh_degree = lp.sh_degree\n",
    "        self.device = gaussians.get_xyz.device\n",
    "        self.compute_cov3D_python = compute_cov3D_python\n",
    "\n",
    "        self.block_dim = lp.block_dim\n",
    "        self.num_cell = lp.block_dim[0] * lp.block_dim[1] * lp.block_dim[2]\n",
    "        self.aabb = lp.aabb\n",
    "        self.scale = scale\n",
    "        self.range = range\n",
    "\n",
    "        self.cell_divider(gaussians)\n",
    "        self.cell_corners = torch.stack(self.cell_corners, dim=0)\n",
    "\n",
    "    def cell_divider(self, gaussians):\n",
    "        with torch.no_grad():\n",
    "            xyz = gaussians.get_xyz\n",
    "            for cell_idx in range(self.num_cell):\n",
    "                cell_mask = block_filtering(cell_idx, xyz, self.aabb, self.block_dim, self.scale)\n",
    "                if self.compute_cov3D_python:\n",
    "                    geometry = gaussians.get_covariance(self.scale)[cell_mask].to(self.device)\n",
    "                else:\n",
    "                    geometry = torch.cat([gaussians.get_scaling[cell_mask],\n",
    "                                          gaussians.get_rotation[cell_mask]], dim=1)\n",
    "\n",
    "                self.feats.append(torch.cat([xyz[cell_mask], \n",
    "                                             gaussians.get_opacity[cell_mask],  \n",
    "                                             gaussians.get_features[cell_mask].reshape(cell_mask.sum(), -1),\n",
    "                                             geometry], dim=1))\n",
    "                xyz_min = xyz[cell_mask].min(dim=0)[0]\n",
    "                xyz_max = xyz[cell_mask].max(dim=0)[0]\n",
    "                corners = torch.tensor([[xyz_min[0], xyz_min[1], xyz_min[2]],\n",
    "                                       [xyz_min[0], xyz_min[1], xyz_max[2]],\n",
    "                                       [xyz_min[0], xyz_max[1], xyz_min[2]],\n",
    "                                       [xyz_min[0], xyz_max[1], xyz_max[2]],\n",
    "                                       [xyz_max[0], xyz_min[1], xyz_min[2]],\n",
    "                                       [xyz_max[0], xyz_min[1], xyz_max[2]],\n",
    "                                       [xyz_max[0], xyz_max[1], xyz_min[2]],\n",
    "                                       [xyz_max[0], xyz_max[1], xyz_max[2]]], device=xyz.device)\n",
    "                self.cell_corners.append(corners)\n",
    "    \n",
    "    def get_feats(self, indices):\n",
    "        out = []\n",
    "        if len(indices) > self.range[0]:\n",
    "            indice_end = min(len(indices), self.range[1])\n",
    "            out = [self.feats[indices[i].item()] for i in range(self.range[0], indice_end)]\n",
    "        return out\n",
    "\n",
    "class BlockedGaussianV3:\n",
    "\n",
    "    gaussians : GaussianModel\n",
    "\n",
    "    def __init__(self, gaussians, lp, range=[0, 1], scale=1.0, compute_cov3D_python=False):\n",
    "        self.cell_corners = []\n",
    "        self.feats = None\n",
    "        self.max_sh_degree = lp.sh_degree\n",
    "        self.device = gaussians.get_xyz.device\n",
    "        self.compute_cov3D_python = compute_cov3D_python\n",
    "        self.cell_ids = torch.zeros(gaussians.get_opacity.shape[0], dtype=torch.long, device=self.device)\n",
    "        self.mask = torch.zeros(gaussians.get_opacity.shape[0], dtype=torch.bool, device=self.device)\n",
    "\n",
    "        self.block_dim = lp.block_dim\n",
    "        self.num_cell = lp.block_dim[0] * lp.block_dim[1] * lp.block_dim[2]\n",
    "        self.aabb = lp.aabb\n",
    "        self.scale = scale\n",
    "        self.range = range\n",
    "\n",
    "        self.cell_divider(gaussians)\n",
    "        self.cell_corners = torch.stack(self.cell_corners, dim=0)\n",
    "\n",
    "    def cell_divider(self, gaussians, n=4):\n",
    "        with torch.no_grad():\n",
    "            if self.compute_cov3D_python:\n",
    "                geometry = gaussians.get_covariance(self.scale).to(self.device)\n",
    "            else:\n",
    "                geometry = torch.cat([gaussians.get_scaling,\n",
    "                                      gaussians.get_rotation], dim=1)\n",
    "            self.feats = torch.cat([gaussians.get_xyz,\n",
    "                                    gaussians.get_opacity,  \n",
    "                                    gaussians.get_features.reshape(geometry.shape[0], -1),\n",
    "                                    geometry], dim=1)\n",
    "\n",
    "            xyz = gaussians.get_xyz\n",
    "            for cell_idx in range(self.num_cell):\n",
    "                cell_mask = block_filtering(cell_idx, self.feats[:, :3], self.aabb, self.block_dim, self.scale)\n",
    "                self.cell_ids[cell_mask] = cell_idx\n",
    "                # MAD to eliminate influence of outsiders\n",
    "                xyz_median = torch.median(xyz[cell_mask], dim=0)[0]\n",
    "                delta_median = torch.median(torch.abs(xyz[cell_mask] - xyz_median), dim=0)[0]\n",
    "                xyz_min = xyz_median - n * delta_median\n",
    "                xyz_min = torch.max(xyz_min, torch.min(xyz[cell_mask], dim=0)[0])\n",
    "                xyz_max = xyz_median + n * delta_median\n",
    "                xyz_max = torch.min(xyz_max, torch.max(xyz[cell_mask], dim=0)[0])\n",
    "                corners = torch.tensor([[xyz_min[0], xyz_min[1], xyz_min[2]],\n",
    "                                       [xyz_min[0], xyz_min[1], xyz_max[2]],\n",
    "                                       [xyz_min[0], xyz_max[1], xyz_min[2]],\n",
    "                                       [xyz_min[0], xyz_max[1], xyz_max[2]],\n",
    "                                       [xyz_max[0], xyz_min[1], xyz_min[2]],\n",
    "                                       [xyz_max[0], xyz_min[1], xyz_max[2]],\n",
    "                                       [xyz_max[0], xyz_max[1], xyz_min[2]],\n",
    "                                       [xyz_max[0], xyz_max[1], xyz_max[2]]], device=xyz.device)\n",
    "                self.cell_corners.append(corners)\n",
    "    \n",
    "    def get_feats(self, indices):\n",
    "        out = torch.tensor([], device=self.device, dtype=self.feats.dtype)\n",
    "        if len(indices) > self.range[0]:\n",
    "            indice_end = min(len(indices), self.range[1])\n",
    "            self.mask = torch.isin(self.cell_ids, indices[self.range[0]:indice_end].to(self.device))\n",
    "            out = self.feats[self.mask]\n",
    "        return out\n",
    "\n",
    "class Camera(nn.Module):\n",
    "    def __init__(self, colmap_id, R, T, FoVx, FoVy, image, gt_alpha_mask,\n",
    "                 image_name, uid,\n",
    "                 trans=np.array([0.0, 0.0, 0.0]), scale=1.0, data_device = \"cuda\"\n",
    "                 ):\n",
    "        super(Camera, self).__init__()\n",
    "\n",
    "        self.uid = uid\n",
    "        self.colmap_id = colmap_id\n",
    "        self.R = R\n",
    "        self.T = T\n",
    "        self.FoVx = FoVx\n",
    "        self.FoVy = FoVy\n",
    "        self.image_name = image_name\n",
    "\n",
    "        try:\n",
    "            self.data_device = torch.device(data_device)\n",
    "        except Exception as e:\n",
    "            print(e)\n",
    "            print(f\"[Warning] Custom device {data_device} failed, fallback to default cuda device\" )\n",
    "            self.data_device = torch.device(\"cuda\")\n",
    "\n",
    "        self.original_image = image.clamp(0.0, 1.0).to(self.data_device)\n",
    "        self.image_width = self.original_image.shape[2]\n",
    "        self.image_height = self.original_image.shape[1]\n",
    "\n",
    "        if gt_alpha_mask is not None:\n",
    "            self.original_image *= gt_alpha_mask.to(self.data_device)\n",
    "        else:\n",
    "            self.original_image *= torch.ones((1, self.image_height, self.image_width), device=self.data_device)\n",
    "\n",
    "        self.zfar = 100.0\n",
    "        self.znear = 0.01\n",
    "\n",
    "        self.trans = trans\n",
    "        self.scale = scale\n",
    "\n",
    "        # world_view_transform = getWorld2View2(R, T, trans, scale)\n",
    "        # if angle_delta is not None:\n",
    "        #     R = world_view_transform[:3, :3]\n",
    "        #     euler = np.array(mat2euler(R)) * 180 / np.pi\n",
    "        #     euler += angle_delta\n",
    "        #     R = euler2mat(*euler * np.pi / 180)\n",
    "        #     world_view_transform[:3, :3] = R\n",
    "\n",
    "        self.world_view_transform = torch.tensor(getWorld2View2(R, T, trans, scale)).transpose(0, 1).cuda()\n",
    "        self.projection_matrix = getProjectionMatrix(znear=self.znear, zfar=self.zfar, fovX=self.FoVx, fovY=self.FoVy).transpose(0,1).cuda()\n",
    "        self.full_proj_transform = (self.world_view_transform.unsqueeze(0).bmm(self.projection_matrix.unsqueeze(0))).squeeze(0)\n",
    "        self.camera_center = self.world_view_transform.inverse()[3, :3]\n",
    "\n",
    "def parse_cfg(cfg):\n",
    "    lp = GroupParams()\n",
    "    op = GroupParams()\n",
    "    pp = GroupParams()\n",
    "\n",
    "    for arg in cfg['model_params'].items():\n",
    "        setattr(lp, arg[0], arg[1])\n",
    "    \n",
    "    for arg in cfg['optim_params'].items():\n",
    "        setattr(op, arg[0], arg[1]) \n",
    "\n",
    "    for arg in cfg['pipeline_params'].items():\n",
    "        setattr(pp, arg[0], arg[1])\n",
    "    \n",
    "    return lp, op, pp\n",
    "\n",
    "def in_frustum(viewpoint_cam, cell_corners, aabb, block_dim):\n",
    "    num_cell = cell_corners.shape[0]\n",
    "    device = cell_corners.device\n",
    "\n",
    "    cell_corners = torch.cat([cell_corners, torch.ones_like(cell_corners[..., [0]])], dim=-1)\n",
    "    cam_center = viewpoint_cam.camera_center\n",
    "    full_proj_transform = viewpoint_cam.full_proj_transform.repeat(num_cell, 1, 1)\n",
    "    viewmatrix = viewpoint_cam.world_view_transform.repeat(num_cell, 1, 1)\n",
    "    cell_corners_screen = cell_corners.bmm(full_proj_transform)\n",
    "    cell_corners_screen = cell_corners_screen / cell_corners_screen[..., [-1]]\n",
    "    cell_corners_screen = cell_corners_screen[..., :-1]\n",
    "\n",
    "    cell_corners_cam = cell_corners.bmm(viewmatrix)\n",
    "    mask = (cell_corners_cam[..., 2] > 0.2)\n",
    "\n",
    "    cell_corners_screen_min = torch.zeros((num_cell, 3), dtype=torch.float32, device=device)\n",
    "    cell_corners_screen_max = torch.zeros((num_cell, 3), dtype=torch.float32, device=device)\n",
    "\n",
    "    for i in range(num_cell):\n",
    "        if mask[i].sum() > 0:\n",
    "            cell_corners_screen_min[i] = cell_corners_screen[i][mask[i]].min(dim=0).values\n",
    "            cell_corners_screen_max[i] = cell_corners_screen[i][mask[i]].max(dim=0).values\n",
    "\n",
    "    box_a = torch.cat([cell_corners_screen_min[:, :2], cell_corners_screen_max[:, :2]], dim=1)\n",
    "    box_b = torch.tensor([[-1, -1, 1, 1]], dtype=torch.float32, device=device)\n",
    "    A = box_a.size(0)\n",
    "    B = box_b.size(0)\n",
    "    max_xy = torch.min(box_a[:, 2:].unsqueeze(1).expand(A, B, 2),\n",
    "                    box_b[:, 2:].unsqueeze(0).expand(A, B, 2))\n",
    "    min_xy = torch.max(box_a[:, :2].unsqueeze(1).expand(A, B, 2),\n",
    "                    box_b[:, :2].unsqueeze(0).expand(A, B, 2))\n",
    "    inter = torch.clamp((max_xy - min_xy), min=0)\n",
    "    mask = (inter[:, 0, 0] * inter[:, 0, 1]) > 0\n",
    "\n",
    "    cam_center_id = which_block(cam_center[None, :], aabb, block_dim)[0]\n",
    "    mask[cam_center_id] = True\n",
    "\n",
    "    return mask\n",
    "\n",
    "def load_gaussians(cfg, config_name, iteration=30_000, load_vq=False, deivce='cuda', source_path='../data/matrix_city/aerial/test/block_all_test'):\n",
    "    \n",
    "    lp, op, pp = parse_cfg(cfg)\n",
    "    setattr(lp, 'config_path', cfg)\n",
    "    lp.source_path = source_path\n",
    "    lp.model_path = os.path.join(\"../output/\", config_name)\n",
    "\n",
    "    modules = __import__('scene')\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        if 'apply_voxelize' in lp.model_config['kwargs'].keys():\n",
    "            lp.model_config['kwargs']['apply_voxelize'] = False\n",
    "        gaussians = getattr(modules, lp.model_config['name'])(lp.sh_degree, device=deivce, **lp.model_config['kwargs'])\n",
    "        scene = LargeScene(lp, gaussians, load_iteration=iteration, load_vq=load_vq, shuffle=False)\n",
    "        print(f'Init {config_name} with {len(gaussians.get_opacity)} points\\n')\n",
    "\n",
    "    return gaussians, scene\n",
    "\n",
    "def loadCam(args, id, cam_info, resolution_scale, angle_delta=None):\n",
    "    orig_w, orig_h = cam_info.image.size\n",
    "\n",
    "    if args.resolution in [1, 2, 4, 8]:\n",
    "        resolution = round(orig_w/(resolution_scale * args.resolution)), round(orig_h/(resolution_scale * args.resolution))\n",
    "    else:  # should be a type that converts to float\n",
    "        if args.resolution == -1:\n",
    "            if orig_w > 1600:\n",
    "                global WARNED\n",
    "                if not WARNED:\n",
    "                    print(\"[ INFO ] Encountered quite large input images (>1.6K pixels width), rescaling to 1.6K.\\n \"\n",
    "                        \"If this is not desired, please explicitly specify '--resolution/-r' as 1\")\n",
    "                    WARNED = True\n",
    "                global_down = orig_w / 1600\n",
    "            else:\n",
    "                global_down = 1\n",
    "        else:\n",
    "            global_down = orig_w / args.resolution\n",
    "\n",
    "        scale = float(global_down) * float(resolution_scale)\n",
    "        resolution = (int(orig_w / scale), int(orig_h / scale))\n",
    "\n",
    "    resized_image_rgb = PILtoTorch(cam_info.image, resolution)\n",
    "\n",
    "    gt_image = resized_image_rgb[:3, ...]\n",
    "    loaded_mask = None\n",
    "\n",
    "    if resized_image_rgb.shape[1] == 4:\n",
    "        loaded_mask = resized_image_rgb[3:4, ...]\n",
    "\n",
    "    if angle_delta is not None:\n",
    "        Rt = np.zeros((4, 4))\n",
    "        Rt[:3, :3] = cam_info.R.transpose()\n",
    "        Rt[:3, 3] = cam_info.T\n",
    "        Rt[3, 3] = 1.0\n",
    "\n",
    "        C2W = np.linalg.inv(Rt)\n",
    "        euler = np.array(mat2euler(C2W[:3, :3])) * 180 / np.pi\n",
    "        euler += angle_delta\n",
    "        C2W[:3, :3] = euler2mat(*euler * np.pi / 180)\n",
    "        Rt = np.linalg.inv(C2W)\n",
    "\n",
    "        R = Rt[:3, :3].transpose()\n",
    "        T = Rt[:3, 3]\n",
    "    else:\n",
    "        R = cam_info.R\n",
    "        T = cam_info.T\n",
    "\n",
    "    return Camera(colmap_id=cam_info.uid, R=R, T=T, \n",
    "                  FoVx=cam_info.FovX, FoVy=cam_info.FovY, \n",
    "                  image=gt_image, gt_alpha_mask=loaded_mask,\n",
    "                  image_name=cam_info.image_name, uid=id, data_device=args.data_device)\n",
    "\n",
    "def loadCamV2(args, id, cam_info, resolution_scale, pitch=None, height=None):\n",
    "    orig_w, orig_h = cam_info.image.size\n",
    "\n",
    "    if args.resolution in [1, 2, 4, 8]:\n",
    "        resolution = round(orig_w/(resolution_scale * args.resolution)), round(orig_h/(resolution_scale * args.resolution))\n",
    "    else:  # should be a type that converts to float\n",
    "        if args.resolution == -1:\n",
    "            if orig_w > 1600:\n",
    "                global WARNED\n",
    "                if not WARNED:\n",
    "                    print(\"[ INFO ] Encountered quite large input images (>1.6K pixels width), rescaling to 1.6K.\\n \"\n",
    "                        \"If this is not desired, please explicitly specify '--resolution/-r' as 1\")\n",
    "                    WARNED = True\n",
    "                global_down = orig_w / 1600\n",
    "            else:\n",
    "                global_down = 1\n",
    "        else:\n",
    "            global_down = orig_w / args.resolution\n",
    "\n",
    "        scale = float(global_down) * float(resolution_scale)\n",
    "        resolution = (int(orig_w / scale), int(orig_h / scale))\n",
    "\n",
    "    resized_image_rgb = PILtoTorch(cam_info.image, resolution)\n",
    "\n",
    "    gt_image = resized_image_rgb[:3, ...]\n",
    "    loaded_mask = None\n",
    "\n",
    "    if resized_image_rgb.shape[1] == 4:\n",
    "        loaded_mask = resized_image_rgb[3:4, ...]\n",
    "\n",
    "    if height is not None and pitch is not None:\n",
    "        Rt = np.zeros((4, 4))\n",
    "        Rt[:3, :3] = cam_info.R.transpose()\n",
    "        Rt[:3, 3] = cam_info.T\n",
    "        Rt[3, 3] = 1.0\n",
    "\n",
    "        C2W = np.linalg.inv(Rt)\n",
    "        euler = np.array(mat2euler(C2W[:3, :3])) * 180 / np.pi\n",
    "        euler[0] = pitch\n",
    "        C2W[:3, :3] = euler2mat(*euler * np.pi / 180)\n",
    "        C2W[2, -1] = height\n",
    "        Rt = np.linalg.inv(C2W)\n",
    "\n",
    "        R = Rt[:3, :3].transpose()\n",
    "        T = Rt[:3, 3]\n",
    "    else:\n",
    "        R = cam_info.R\n",
    "        T = cam_info.T\n",
    "\n",
    "    return Camera(colmap_id=cam_info.uid, R=R, T=T, \n",
    "                  FoVx=cam_info.FovX, FoVy=cam_info.FovY, \n",
    "                  image=gt_image, gt_alpha_mask=loaded_mask,\n",
    "                  image_name=cam_info.image_name, uid=id, data_device=args.data_device)\n",
    "\n",
    "def render_lod(viewpoint_cam, lod_list : list, pipe, bg_color : torch.Tensor, scaling_modifier = 1.0, override_color = None):\n",
    "        \n",
    "        # sort cells by distance to camera\n",
    "\n",
    "        in_frustum_mask = in_frustum(viewpoint_cam, lod_list[-1].cell_corners, lod_list[-1].aabb, lod_list[-1].block_dim)\n",
    "        in_frustum_indices = in_frustum_mask.nonzero().squeeze(0)\n",
    "        cam_center = viewpoint_cam.camera_center\n",
    "        distance3D = torch.norm(lod_list[-1].cell_corners[in_frustum_mask, :, :2] - cam_center[:2], dim=2).min(dim=1).values\n",
    "        in_frustum_indices = in_frustum_indices[torch.sort(distance3D)[1]]\n",
    "        \n",
    "        # used for BlockedGaussianV3\n",
    "        out_list = []\n",
    "        main_device = lod_list[-1].feats.device\n",
    "        max_sh_degree = lod_list[-1].max_sh_degree\n",
    "        feat_end_dim = 3 * (max_sh_degree + 1) ** 2 + 4\n",
    "        \n",
    "        for lod_gs in lod_list:\n",
    "            out_i = lod_gs.get_feats(in_frustum_indices)\n",
    "            if out_i.shape[0] == 0:\n",
    "                continue\n",
    "            if out_i.device != main_device:\n",
    "                out_i = torch.cat([out_i[:, :3].to(main_device), out_i[:, 3:].half().to(main_device)], dim=1)\n",
    "            out_list.append(out_i)\n",
    "\n",
    "        feats = torch.cat(out_list, dim=0)\n",
    "        # feats = lod_list[1].feats\n",
    "\n",
    "        means3D = feats[:, :3].float()\n",
    "        screenspace_points = torch.zeros_like(means3D, dtype=means3D.dtype, requires_grad=True, device=\"cuda\") + 0\n",
    "        means2D = screenspace_points\n",
    "        opacity = feats[:, 3].float()\n",
    "        # If precomputed 3d covariance is provided, use it. If not, then it will be computed from\n",
    "        # scaling / rotation by the rasterizer.\n",
    "        scales = None\n",
    "        rotations = None\n",
    "        cov3D_precomp = None\n",
    "        if pipe.compute_cov3D_python:\n",
    "            cov3D_precomp = feats[:, feat_end_dim:].float()\n",
    "        else:\n",
    "            scales = feats[:, feat_end_dim:feat_end_dim+3].float()\n",
    "            rotations = feats[:, (feat_end_dim+3):].float()\n",
    "            \n",
    "        # If precomputed colors are provided, use them. Otherwise, if it is desired to precompute colors\n",
    "        # from SHs in Python, do it. If not, then SH -> RGB conversion will be done by rasterizer.\n",
    "        shs = None\n",
    "        colors_precomp = None\n",
    "        if override_color is None:\n",
    "            features = feats[:, 4:feat_end_dim].reshape(-1, (max_sh_degree+1)**2, 3).float()\n",
    "            if pipe.convert_SHs_python:\n",
    "                shs_view = features.transpose(1, 2).view(-1, 3, (max_sh_degree+1)**2)\n",
    "                dir_pp = (means3D - viewpoint_cam.camera_center.repeat(features.shape[0], 1))\n",
    "                dir_pp_normalized = dir_pp/dir_pp.norm(dim=1, keepdim=True)\n",
    "                sh2rgb = eval_sh(max_sh_degree, shs_view, dir_pp_normalized)\n",
    "                colors_precomp = torch.clamp_min(sh2rgb + 0.5, 0.0)\n",
    "            else:\n",
    "                shs = features\n",
    "        else:\n",
    "            colors_precomp = override_color  # check if requires masking\n",
    "        \n",
    "        # Set up rasterization configuration\n",
    "        tanfovx = math.tan(viewpoint_cam.FoVx * 0.5)\n",
    "        tanfovy = math.tan(viewpoint_cam.FoVy * 0.5)\n",
    "\n",
    "        raster_settings = GaussianRasterizationSettings(\n",
    "            image_height=int(viewpoint_cam.image_height),\n",
    "            image_width=int(viewpoint_cam.image_width),\n",
    "            tanfovx=tanfovx,\n",
    "            tanfovy=tanfovy,\n",
    "            bg=bg_color,\n",
    "            scale_modifier=scaling_modifier,\n",
    "            viewmatrix=viewpoint_cam.world_view_transform,\n",
    "            projmatrix=viewpoint_cam.full_proj_transform,\n",
    "            sh_degree=max_sh_degree,\n",
    "            campos=viewpoint_cam.camera_center, \n",
    "            prefiltered=False,\n",
    "            debug=pipe.debug\n",
    "        )\n",
    "\n",
    "        rasterizer = GaussianRasterizer(raster_settings=raster_settings)\n",
    "\n",
    "        rendered_image, radii = rasterizer(\n",
    "            means3D = means3D,\n",
    "            means2D = means2D,\n",
    "            shs = shs,\n",
    "            colors_precomp = colors_precomp,\n",
    "            opacities = opacity,\n",
    "            scales = scales,\n",
    "            rotations = rotations,\n",
    "            cov3D_precomp = cov3D_precomp)\n",
    "        \n",
    "        # Those Gaussians that were frustum culled or had a radius of 0 were not visible.\n",
    "        # They will be excluded from value updates used in the splitting criteria.\n",
    "        return {\"render\": rendered_image,\n",
    "                \"viewspace_points\": screenspace_points,\n",
    "                \"visibility_filter\" : radii > 0,\n",
    "                \"radii\": radii}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Render the whole dataset without LoD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import torchvision\n",
    "from os import makedirs\n",
    "from torch.profiler import profile, record_function, ProfilerActivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading camera 741/741\n",
      "Init block_mc_aerial_block_all_lr_c36_loss_5_50_lr64_vq with 11782801 points\n",
      "\n"
     ]
    }
   ],
   "source": [
    "iteration = 30_000\n",
    "\n",
    "config = '../config/block_mc_aerial_block_all_lr_c36_loss_5_50_lr64_vq.yaml'\n",
    "config_name = os.path.splitext(os.path.basename(config))[0]\n",
    "with open(config) as f:\n",
    "    cfg = yaml.load(f, Loader=yaml.FullLoader)\n",
    "    # cfg['model_params']['model_config']['name'] = \"GaussianModel\"\n",
    "gaussians, scene = load_gaussians(cfg, config_name, iteration=None, load_vq=True, deivce='cuda')\n",
    "\n",
    "lp, op, pp = parse_cfg(cfg)\n",
    "setattr(lp, 'config_path', cfg)\n",
    "lp.source_path = '../data/matrix_city/aerial/test/block_all_test'\n",
    "scaling_modifier = 1.0\n",
    "override_color = None\n",
    "angle_delta = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_path = os.path.join(\"../output/\", \"block_mc_aerial_block_all_lr_c36_loss_5\")\n",
    "name = \"block_all_test\"\n",
    "\n",
    "# indexing is the most time_consuming part\n",
    "# angle_delta = np.array([0, 0, 0])\n",
    "# render_path = os.path.join(model_path, name, \"ours_{}\".format(iteration), \"renders\")\n",
    "# gts_path = os.path.join(model_path, name, \"ours_{}\".format(iteration), \"gt\")\n",
    "\n",
    "pitch, height = -180, 20\n",
    "# render_path = os.path.join(model_path, name, \"ours_{}_pitch{}_height{}\".format(iteration, pitch, height), \"renders\")\n",
    "# gts_path = os.path.join(model_path, name, \"ours_{}_pitch{}_height{}\".format(iteration, pitch, height), \"gt\")\n",
    "\n",
    "# makedirs(render_path, exist_ok=True)\n",
    "# makedirs(gts_path, exist_ok=True)\n",
    "pynvml.nvmlInit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 741/741 [02:22<00:00,  5.20it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average FPS: 57.2602\n",
      "Min FPS: 1.1795\n",
      "Average Memory: 7312.5326 M\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    \n",
    "    views = scene.getTrainCameras()  # getTrainCameras, getTestCameras\n",
    "    bg_color = [1,1,1] if lp.white_background else [0, 0, 0]\n",
    "    background = torch.tensor(bg_color, dtype=torch.float32, device=\"cuda\")\n",
    "    avg_render_time = 0\n",
    "    max_render_time = 0\n",
    "    avg_memory = 0\n",
    "    visible_pts_cnt_single = []\n",
    "\n",
    "    for idx in tqdm(range(len(views))):\n",
    "        \n",
    "    # idx = 518\n",
    "\n",
    "        # viewpoint_cam = loadCam(lp, idx, views[idx], 1.0, angle_delta)\n",
    "        viewpoint_cam = loadCamV2(lp, idx, views[idx], 1.0, pitch, height)\n",
    "\n",
    "        torch.cuda.empty_cache()\n",
    "\n",
    "        start = time.time()\n",
    "        render_pkg = render(viewpoint_cam, gaussians, pp, background)\n",
    "        end = time.time()\n",
    "\n",
    "        handle = pynvml.nvmlDeviceGetHandleByIndex(5)\n",
    "        memory_info = pynvml.nvmlDeviceGetMemoryInfo(handle)\n",
    "        memory_used = memory_info.used / 1024 / 1024\n",
    "        visible_pts_cnt_single.append(render_pkg['visibility_filter'].sum().item())\n",
    "\n",
    "        # torch.cuda.empty_cache()\n",
    "\n",
    "    # image = render_pkg[\"render\"]\n",
    "\n",
    "    # image = image.cpu().numpy().transpose(1,2,0)\n",
    "\n",
    "    # show render results\n",
    "    # plt.figure(figsize=(8, 8))\n",
    "    # plt.imshow(image)\n",
    "    # plt.title(f\"Idx {idx}, Rendered {render_pkg['visibility_filter'].sum().item()} / {render_pkg['radii'].shape[0]} points\")\n",
    "    # plt.axis(False)\n",
    "    # plt.tight_layout()\n",
    "    # plt.show()\n",
    "        avg_memory += memory_used\n",
    "        avg_render_time += end-start\n",
    "        max_render_time = max(max_render_time, end-start)\n",
    "\n",
    "        # torchvision.utils.save_image(render_pkg[\"render\"], os.path.join(render_path, '{0:05d}'.format(idx) + \".png\"))\n",
    "        # torchvision.utils.save_image(viewpoint_cam.original_image[0:3, :, :], os.path.join(gts_path, '{0:05d}'.format(idx) + \".png\"))\n",
    "\n",
    "    print(f'Average FPS: {len(views)/avg_render_time:.4f}')\n",
    "    print(f'Min FPS: {1/max_render_time:.4f}')\n",
    "    print(f'Average Memory: {avg_memory/len(views):.4f} M')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAHACAYAAABwEmgAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABC6ElEQVR4nO3deVhU5f//8ReIS6BiKmrmrjFoLOISqShJGuaWuJapaWRqomb6yz0zNNRccl/STCtTyqVQtLI+mpZbfSgzyfqUAmoh4oKCCcj8/vBivo0sOshwQJ+P6+K6mPucM+d9D8PMa+5znzMOZrPZLAAAAIM4Gl0AAAC4txFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUZgmIMHD8pkMmnnzp1Gl2KTzZs3y2Qy6dSpU3bf1/jx4xUYGGi5ferUKZlMJq1evdru+5akRYsWyWQyFcq+8isjI0OzZ89WQECAPDw89NJLLxXI/fbv31/9+/e3aZucnhuBgYEaMmTILbfN+n84ePCgzbXe7bKe95s3bza6FNgJYeQul/Xi6OXlpYSEhGzL+/fvr86dOxtQWdGQ9QaQ9ePp6amWLVuqf//+Wr58uc6fP18g+7l69aoWLVpUJN9oinJtt2PTpk1avXq1goKCNHPmTA0cODDbOklJSWrUqJHGjh2b6/1cuXJF3t7eCg0NtWO1he/KlStatmyZunfvrqZNm8rT01Nt27bVyy+/rN27dxtdHiBJcjK6ABSOtLQ0rVy5UlOmTDG6lCKpf//+8vLyUmZmps6fP6/o6GgtWrRIa9as0dtvv60WLVpY1n3qqafUqVMnlSpV6rbv/+rVq1q8eLFCQ0Pl5+d329uFhYXJ3l8flVdtw4YN04svvmjX/d+pAwcOqGrVqpo4cWKu61SqVEktW7bUV199patXr+q+++7Lts6XX36pa9euqWvXrpKUr9Gn/Dw37Ck2NlYhISE6c+aM2rVrp27dusnZ2Vl///239uzZoyFDhmjWrFnq1q2b0aXm6cEHH9SRI0fk5MRb1t2Kv+w9omHDhoqIiNCLL76oqlWrGl1OoUpNTZWzs3Oe6zRr1kwdOnSwavv111/1/PPPa+TIkdq+fbuqVKkiSSpRooRKlChht3ql/6u5ZMmSdt3PrTg5ORX5N4CkpCSVL1/+lut16dJFe/fu1ddff61OnTplW75t2zaVK1dOjz32mCTlK1AUxnPjdmVkZCg0NFRJSUl6//331bRpU6vloaGh2rdvn65fv25QhbfPwcFBpUuXNroM2BGHae4RQ4YMUWZmpt55550818vr2KzJZNKiRYsst7PmE5w4cUJjx45V06ZN9eijj+rtt9+W2WzWX3/9pWHDhqlJkyZq1aqV3n333Rz3mZmZqXnz5qlVq1Zq3Lixhg4dqr/++ivbej/99JNCQkLUtGlT+fj4qF+/fvrhhx+s1smq6X//+5/GjBmj5s2bq2/fvrfzEGXj4eGhiRMnKjk5WR9++KGlPad5AT///LNCQkLk5+cnb29vBQYGasKECZJuPKZZIyuLFy+2HBLKeizHjx8vX19fxcXFafDgwfL19bUcTrh5zsi/vffee2rbtq28vb3Vr18//fbbb1bLc5vz8O/7vFVtOc0ZycjI0JIlS9SuXTt5enoqMDBQ8+bNU1pamtV6WXMlvv/+e/Xs2VNeXl56/PHHtXXr1jwe9f+TmpqqmTNnKiAgQJ6engoKCtLq1astI0VZz9WDBw/q999/t9Se2+Gm9u3by9nZWZGRkdmWJSUlaf/+/QoKCrKEkJwev/fff1+dOnWSj4+Pmjdvru7du1vdX17zifbt26ennnpKXl5e6tixo7744ovbehxu53mfk507d+q3337TsGHDsgWRLP7+/goICLDcvnjxombNmqUuXbrI19dXTZo00QsvvKBff/3Varvc+pnTvJeTJ09qxIgRatWqlby8vNSmTRuNHj1aly9ftqzz7bff6plnnlGzZs3k6+uroKAgzZs3z7I8p9elX3/9VePHj9fjjz8uLy8vtWrVShMmTNCFCxesasp6DsfGxmr8+PFq1qyZmjZtqgkTJujq1atW696qDthP0f7IgwJTo0YNPfXUU4qIiNDgwYMLdHRk9OjRql+/vsaMGaM9e/Zo2bJlqlChgjZs2KBHH31UY8eOVWRkpGbNmiUvLy81b97cavtly5bJwcFBgwcPVlJSktauXauBAwfq008/VZkyZSRJ+/fv1+DBg+Xp6anQ0FA5ODho8+bNeu6557R+/Xp5e3tb3eeoUaNUu3ZtjR49+o4OcwQFBWnSpEnat2+fRo8eneM6SUlJCgkJ0f33368XX3xR5cuX16lTp/Tll19KkipWrKjXX39dr7/+utq3b6/27dtLktWbfEZGhuUNZ9y4cZZ+52br1q1KSUlR3759de3aNb3//vt67rnnFBkZqcqVK992/26ntptNnjxZW7ZsUVBQkAYNGqQjR45oxYoV+uOPP7RkyRKrdWNjYzVq1Cj17NlTwcHB2rRpk8aPH6+HH35YDz30UK77MJvNGjZsmA4ePKiePXuqYcOG2rt3r2bPnq2EhARNnDhRFStW1OzZs7V8+XKlpqbqlVdekSTVr18/x/t0dnZWYGCgPv/8c128eFEVKlSwLIuKitL169fVpUuXXGuKiIjQ9OnTFRQUpAEDBujatWs6fvy4fvrppzy3k268IY8ePVpPP/205XEYNWqUVq1apVatWuW6na3P+3/7z3/+I+nGoaPbFR8fr127dqlDhw6qUaOGzp07p40bN6pfv37avn27za8baWlpCgkJUVpamvr166fKlSsrISFBu3fvVnJyssqVK6fff/9dQ4YMkclk0siRI1WqVCnFxsbqv//9b573/d133yk+Pl7du3eXm5ubfv/9d0VEROh///ufIiIi5ODgYLX+yy+/rBo1auiVV17RsWPH9PHHH6tixYr6f//v/0lSvutAATHjrrZp0yazu7u7+ciRI+a4uDhzo0aNzGFhYZbl/fr1M3fq1MlyOz4+3uzu7m7etGlTtvtyd3c3L1y40HJ74cKFZnd3d/OUKVMsbRkZGeY2bdqYTSaTecWKFZb2S5cumb29vc3jxo2ztB04cMDs7u5ubt26tfny5cuW9qioKLO7u7t57dq1ZrPZbM7MzDQ/8cQT5ueff96cmZlpWe/q1avmwMBA86BBg7LV9Morr9zW45NVw44dO3Jdp2vXrubmzZtbbmc9pvHx8Waz2Wz+8ssvLY9xbpKSkrI9flnGjRtndnd3N8+ZMyfHZW3btrXczvr7eHt7m//++29L+08//WR2d3c3v/nmm5a2fv36mfv163fL+8yrtqzHM0tMTIzZ3d3dPGnSJKv1Zs6caXZ3dzfv37/f0ta2bVuzu7u7+fDhw1b78vT0NM+cOTPbvv4t6zFdunSpVfuIESPMJpPJHBsba9XPfz+H87J7926zu7u7ecOGDVbtvXv3Nrdu3dp8/fp1q/v99+M3bNiwW+7n5ueG2fx/j8Pnn39uabt8+bK5VatW5m7dulnasp6LBw4cMJvNtj3vc9KtWzdzs2bNsrWnpKSYk5KSLD///t+7du2a1WNgNt94znl6epoXL16cZz9z6sOxY8du+f+1Zs0as7u7uzkpKSnXdXJ6Xbp69Wq29bZt25btOZf1HJ4wYYLVusOHDzc/8sgjNtUB++EwzT2kZs2a6tq1qyIiInT27NkCu9+ePXtafi9RooQ8PT1lNput2suXL6+6desqPj4+2/bdunVT2bJlLbc7dOggNzc37dmzR5IUExOjkydPqkuXLrpw4YLOnz+v8+fPKzU1VS1atNDhw4eVmZlpdZ9PP/10gfXP2dlZKSkpuS4vV66cJGn37t1KT0/P936eeeaZ2163Xbt2Vp9Svb295ePjY3nM7CXr/gcNGmTV/vzzz1stz9KgQQM1a9bMcrtixYq5Pg/+7ZtvvlGJEiWyHSZ5/vnnZTab9c033+Sr/latWqlixYratm2bpS0+Pl4//vijOnXqJEfH3F8Sy5cvr7///ltHjhyxeb9VqlSxjDpJUtmyZdWtWzcdO3ZMiYmJOW6Tn+f9v125ciXHuVLz589XixYtLD9jxoyxLCtVqpTlMbh+/bouXLggZ2dn1a1bV8eOHbO531n/1/v27ct2SCRL1nyfr776Ks/+3Ozfo4fXrl3T+fPn5ePjI0n65Zdfsq1/82tCs2bNdPHiRV25cuWO6kDB4DDNPeall17SZ599ppUrV2ry5MkFcp/Vq1e3ul2uXDmVLl1aFStWzNZ+8eLFbNvXrl3b6raDg4Nq166t06dPS7oxxC1J48aNy7WGy5cvy9XV1XK7Ro0atnQhT6mpqXJxccl1+SOPPKKgoCAtXrxY7733nh555BG1a9dOXbp0ue1JkE5OTqpWrdpt13TzYyZJderU0Y4dO277PvLj9OnTcnR0VK1ataza3dzcVL58ecvfLMsDDzyQ7T5cXV116dKlW+6nSpUqViFV+r9DMDfv53Y5OTmpY8eOWr9+vRISElS1alVLMMk6iyY3gwcP1nfffadevXqpdu3aatWqlTp37pzrfIx/q127drbDBnXq1LH0xc3NLds2+Xne/5uLi0uO/299+/ZV27ZtJclyiCJLZmam1q1bp/Xr1+vUqVNWk1v/fVjrdtWsWVODBg3SmjVrFBkZqWbNmikwMFBdu3a1hPiOHTvq448/1uTJkzV37ly1aNFC7du3V4cOHfIMhxcvXtTixYsVFRWlpKQkq2X/no+S5ebXqazwcenSJZUtWzbfdaBgEEbuMf8eHcnplM2bXzCz5DXjPqd/1NzOKDDnY/5G1javvvqqGjZsmOM6N38CLKiZ9+np6Tp58mSe8xscHBy0cOFC/fjjj/rPf/6jvXv3auLEiVqzZo02btyYZ5DJ8u9PpPZWEGdP5PY8uVlRObPk37p27aoPPvhA27ZtU0hIiLZv364GDRrk+tzKUr9+fe3cuVO7d+/W3r179cUXX2j9+vUaPny4Ro4cWeB15ud5/2/16tVTTEyMJXRlqVu3rurWrSsp+//J8uXLtWDBAvXo0UOjRo2Sq6urHB0d9eabb1r97+b2989pRGH8+PEKDg7WV199pW+//VbTp0/XihUrFBERoWrVqqlMmTL68MMPdfDgQctjGxUVpY0bN+rdd9/N9Tn08ssvKzo6WiEhIWrYsKGcnZ2VmZmpF154IcfXmdz+v7LWzW8dKBiEkXvQsGHD9Nlnn+V4Zk3Wp6zk5GSr9jNnztitntjYWKvbZrNZsbGxlkmUNWvWlHRjyLdly5Z2qyMnn3/+uf755x/5+/vfct3GjRurcePGGj16tCIjIzV27FhFRUWpV69et/3mfbtufsykG5+kH3zwQcttV1fXHA+H3Py3tKW2Bx98UJmZmYqNjbWaKHru3DklJydb7f9OPPjgg9q/f7+uXLliNTry559/Wpbnl4+Pj2rVqqVt27apVatW+v3333OdnHwzZ2dndezYUR07dlRaWppGjBih5cuXa8iQIXkG4NjYWJnNZqvHOmvkI7e+3Onz/rHHHtP27dv12WefafDgwbe1zeeffy4/Pz+9+eabVu3Jycm6//77LbezRhVuHoHIbcQq60ynl156Sf/973/1zDPP6KOPPrI87o6OjpbDRhMmTNDy5cs1f/58HTx4MMe+X7p0Sfv379eIESOsLlKX9Zjml611oOAw9nQPqlWrlrp27aqNGzdmO15dtmxZ3X///fr++++t2tevX2+3erZu3Wo5bivdOCUxMTFRbdq0kSR5enqqVq1aevfdd3Ocu1FQV0m92a+//qo333xTrq6uevbZZ3Nd79KlS9k+iWV9ks063TXrIls3h7z82rVrl9UVdY8cOaKffvrJ8phJN97M/vzzT6vH59dff812doAttWWdBrp27Vqr9jVr1lgtv1Nt2rTR9evXrU6plm6czuzg4GDVz/zo0qWLjh07poULF8rBweG2rkJ88ymjpUqVUv369WU2m285V+js2bOWs6ukG/M5tm7dqoYNG+Z4iEa68+f9k08+qQYNGmjp0qX68ccfc1zn5udtiRIlsrXt2LEj29Wbsw7THT582NJ2/fp1RUREWK135coVZWRkWLW5u7vL0dHR8r+R06Gkm/9/bpbbKMXNz0tb5KcOFBxGRu5RQ4cO1aeffqoTJ05kOwTRq1cvrVy5UpMmTZKnp6e+//57nThxwm61uLq6qm/fvurevbvl1N7atWurd+/ekm58Wpk+fboGDx6szp07q3v37qpataoSEhJ08OBBlS1bVsuXL7+jGr7//ntdu3ZNmZmZunjxov773//q66+/VtmyZbV48eJc3zAkacuWLfroo4/Url071apVSykpKYqIiFDZsmUtb5plypRRgwYNtGPHDtWpU0cVKlTQQw89JHd393zVW6tWLT3zzDN65plnlJaWpnXr1qlChQp64YUXLOv07NlT7733nkJCQtSzZ08lJSVpw4YNatCggdWbmy21eXh4KDg4WBs3blRycrKaN2+un3/+WVu2bFG7du306KOP5qs/NwsMDJSfn5/mz5+v06dPy2Qy6dtvv9VXX32l5557LtucFVt17dpVS5Ys0VdffaUmTZrc1hyjkJAQVa5cWU2aNFGlSpX0559/6oMPPlBAQEC2uS03q1OnjiZNmqSff/5ZlSpV0qZNm5SUlKTw8PBct7nT533JkiW1ePFihYSEqG/fvmrfvr2aNWum++67TwkJCfr666915swZqwD52GOPacmSJZowYYJ8fX3122+/KTIy0jJKk+Whhx5S48aNNW/ePF26dEmurq6KiorKFjwOHDigN954Qx06dFCdOnV0/fp1ffrppypRooSCgoIkSUuWLNH333+vgIAAPfjgg0pKStL69etVrVq1XOfjlC1bVs2bN9eqVauUnp6uqlWr6ttvv72j74vKTx0oOISRe1Tt2rXVtWtXbdmyJduy4cOH6/z58/r888+1Y8cOtWnTRqtWrbK6JHpBGjp0qI4fP66VK1cqJSVFLVq00NSpU60u2e3n56eNGzdq6dKl+uCDD5Samio3Nzd5e3urT58+d1zD+++/L+nGC3i5cuVUv359jRgxQr179842EfdmjzzyiH7++WdFRUXp3LlzKleunLy9vTVnzhyrF/Hp06crLCxM4eHhSk9PV2hoaL7DSLdu3eTo6Ki1a9cqKSlJ3t7emjJliuUqsdKNOQ6zZs3SwoULFR4ergYNGmj27Nnatm2bDh06ZHV/ttQ2ffp01ahRQ1u2bNGuXbtUuXJlDRkypEC/08XR0VHLli3TwoULFRUVpc2bN+vBBx/Uq6++ajlz507UqVNHXl5e+vnnn295jZAsffr0UWRkpNasWaPU1FRVq1ZN/fv3v60v5qtTp46mTJmi2bNn68SJE6pRo4bmz5+v1q1b57ndnT7v69atq08//VTr1q3Trl279M033yg9PV2VK1e2fA9P1mRW6cb/4tWrVxUZGamoqCg1atRIK1as0Ny5c7Pd95w5c/Taa69p5cqVKl++vHr27Ck/Pz+rM61MJpP8/f31n//8RwkJCbrvvvtkMpn0zjvvqHHjxpJuBM/Tp09r06ZNunDhgu6//3498sgjGjFihGWSa07mzp2rsLAwrV+/XmazWa1atdI777xzy8c0N/mtAwXDwZyfGYUAAAAFhDkjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMVeRP7c3IyNClS5dUunRpvh8AAIBiIjMzU9euXZOrq6ucnPKOG0U+jFy6dOmOL/ELAACMUadOHVWqVCnPdYp8GMn6voc6depYXQQLAAAUXVevXtXJkydv64tLi3wYyTo0c9999+X5DZUAAKDouZ0pFkzCAAAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChivx30wAAir7aCwp3f7GjCnd/sC/CCADgnjB+/HglJydr6dKl+dp2y5YtkiQnJye5urrKZDKpU6dO6t69+219GRxyRxgBAOA2tG7dWuHh4crMzNS5c+e0d+9ezZgxQ59//rmWLVsmJyfeUvOLKAcAuOcdOnRIPXv2lKenp/z9/TVnzhxlZGRYrVOqVCm5ubmpatWqevjhhzV06FAtXbpU33zzjWXUBPlDjDNQbsdYORYKAIUnISFBL774ooKDgzVr1iydOHFCkydPVunSpTVixIg8t23RooU8PDz0xRdfqFevXoVU8d2HMAIAuKetX79e1apV02uvvSYHBwfVr19fCQkJmjNnjoYPH37L+SD16tXT8ePHC6nau5PNh2kOHz6soUOHyt/fXyaTSbt27bJabjKZcvxZtWqVZZ3AwMBsy1euXHnnvQEAwEZ//PGHfH195eDgYGlr2rSpUlNT9ffff99ye7PZbLUtbGfzyEhqaqpMJpN69Oih0NDQbMv37dtndfubb77RpEmTFBQUZNU+cuRI9e7d23LbxcXF1lIAADDcH3/8oRo1ahhdRrFmcxgJCAhQQEBArsvd3Nysbn/11Vfy8/NTzZo1rdpdXFyyrQsAQGGrX7++Pv/8c6sRjh9++EEuLi6qVq1antvu379fv/32mwYOHFgIld697Dpn5Ny5c9qzZ49mzpyZbdk777yjZcuW6YEHHlDnzp01cOBATosCANjV5cuXFRMTY9XWu3dvrV27VmFhYXr22Wd14sQJLVq0SIMGDbKaL5KWlqbExESrU3tXrFihtm3bqlu3boXck7uLXd/9t2zZIhcXFz3xxBNW7f3791ejRo3k6uqq6OhozZs3T4mJiZowYYI9ywEA2ElxOQvw0KFD2YJDz549tXLlSs2ePVsRERGqUKGCevbsqWHDhlmtt3fvXvn7+8vJyUnly5eXh4eHJk+erODgYC56docczGazOb8bm0wmLVmyRO3atctxeYcOHdSqVStNmTIlz/v55JNPNHXqVEVHR6tUqVJWy1JTUxUTE6OGDRvK2dk5v6UWSZzaCwC4W9ny/m23KPf999/rxIkTt3XetY+PjzIyMnTq1Cl7lQMAAIoou4WRTz75RA8//LA8PDxuuW5MTIwcHR1VqVIle5UDAACKKJvnjKSkpCguLs5y+9SpU4qJiZGrq6uqV68uSbpy5Yp27typcePGZds+OjpaP/30kx599FG5uLgoOjpa4eHh6tq1q1xdXe+gKwAAoDiyOYwcPXpUAwYMsNwODw+XJAUHB1vOmtm+fbvMZrM6d+6cbftSpUopKipKixcvVlpammrUqKGBAwdq0KBB+e0DAAAoxu5oAmthYAIrAADFT5GYwAoAAHA7CCMAAMBQhBEAAGAowggAADAUXwYDALhzrwcX8v62FPhd3uqq4vmxaNEi7dq1S59++mmB3efdiDACALjrnT9/XgsWLNCePXt07tw5ubq6ysPDQy+99JKaNm0qSdq3b1+RvN7V+PHjlZycrKVLl+Z7+y1bboQ3Jycnubq6ymQyqVOnTurevbvle3UuXryoRYsWad++ffrrr79UsWJFtWvXTqNGjVK5cuUKrD85IYwAAO56I0aMUHp6umbOnKmaNWsqKSlJ+/fv18WLFy3ruLm5GVegnbVu3Vrh4eFW3zg8Y8YMff7551q2bJmcnJx09uxZnT17VuPGjVODBg10+vRpvf766zp79qwWLlxo1/qYMwIAuKslJyfr+++/19ixY/Xoo4/qwQcflLe3t4YMGaLHH3/csp7JZNKuXbsk3bi6uMlk0hdffKH+/fvLx8dHXbt2VXR0tNV9R0REKCAgQD4+Pho+fLjWrFmjZs2a5VnPxx9/rCeffFJeXl7q0KGDPvzwwzvq36FDh9SzZ095enrK399fc+bMUUZGhtU6pUqVkpubm6pWraqHH35YQ4cO1dKlS/XNN99YRk3c3d21aNEiBQYGqlatWmrRooVefvllff3119nur6ARRgAAdzVnZ2c5Oztr165dSktLs2nb+fPnKyQkRFu3blWdOnU0ZswYyxvzDz/8oKlTp2rAgAHaunWrWrZsqeXLl+d5f5999pkWLFig0aNHKyoqSq+88ooWLlxoCQS2SkhI0IsvvigvLy99+umnev311/XJJ59o2bJlt9y2RYsW8vDw0BdffJHrOleuXFHZsmXl5GTfAykcpgEA3NWcnJw0c+ZMTZkyRRs2bFCjRo30yCOPqGPHjrf8Mtfnn39ejz32mCRp5MiR6tSpk2JjY1W/fn198MEHatOmjUJCQiRJdevWVXR0tHbv3p3r/S1atEjjx4/XE088IUmqWbOm/ve//2njxo0KDrZ9EvD69etVrVo1vfbaa3JwcFD9+vWVkJCgOXPmaPjw4Zb5ILmpV6+ejh8/nuOy8+fPa+nSperTp4/NddmKkREAwF0vKChIe/fu1bJly9S6dWsdOnRI3bt31+bNm/PczmQyWX7PmlNy/vx5SdKJEyfk5eVltb63t3eu95Wamqq4uDhNmjRJvr6+lp9ly5ZZfQGtLf744w/5+vrKwcHB0ta0aVOlpqbq77//vuX2ZrPZatssV65c0ZAhQ1S/fn2FhobmqzZbMDICALgnlC5dWq1atVKrVq00fPhwTZo0SYsWLVL37t1z3aZkyZKW37PetDMzM/O1/9TUVElSWFiYfHx8rJbdagTDXv744w/VqFHDqu3KlSt64YUX5OLioiVLllg9BvbCyAgA4J7UoEEDS0DIj7p16+ro0aNWbT///HOu61euXFlVqlRRfHy8ateubfVTs2bNfNVQv359RUdH69/fefvDDz/IxcVF1apVy3Pb/fv367fffrMcMpJuBJGQkBCVLFlSy5YtU+nSpfNVl60YGQEA3NUuXLigUaNGqUePHjKZTHJxcdHRo0e1atUqq7NpbNWvXz/169dPa9asUdu2bXXgwAF98803OR72yDJy5EhNnz5d5cqVU+vWrZWWlqajR48qOTlZgwYNynW7y5cvKyYmxqqtQoUK6tu3r9auXauwsDA9++yzOnHihBYtWqRBgwZZjbakpaUpMTHR6tTeFStWqG3bturWrZukG0Hk+eef19WrV/XWW2/pypUrunLliiSpYsWKKlGiRL4fq1shjAAA7pwdrohaUFxcXOTj46O1a9cqLi5OGRkZqlatmnr16qWhQ4fm+36bNm2qadOmafHixXr77bfl7++vgQMH5nmqbq9evVSmTBmtXr1as2fPlrOzs9zd3fXcc8/lua9Dhw5ZQkOWnj17asaMGVq5cqVmz56tiIgIVahQQT179tSwYcOs1t27d6/8/f3l5OSk8uXLy8PDQ5MnT1ZwcLAltPzyyy/66aefJEnt27e32v6rr77KdjinIDmY/z22UwSlpqYqJiZGDRs2lLOzs9HlFKjaC3Jujx1VuHUAAArG5MmT9eeff2r9+vVGl2I4W96/mTMCAEA+rV69Wr/++qtiY2P1/vvva+vWrfk6Rfdex2EaAADy6ciRI1q1apVSUlJUs2ZNTZo0Sb169TK6rGKHMAIAQD4tWJDL8XbYhMM0AADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGMrmMHL48GENHTpU/v7+MplM2rVrl9Xy8ePHy2QyWf2EhIRYrXPx4kWNGTNGTZo0UbNmzTRx4kSlpKTcWU8AAECx5GTrBqmpqTKZTOrRo4dCQ0NzXKd169YKDw+33C5VqpTV8rFjxyoxMVFr1qxRenq6Jk6cqNdee01z5861tRwAAFDM2RxGAgICFBAQkOc6pUqVkpubW47L/vjjD+3du1effPKJvLy8JEmTJ0/Wiy++qFdffVVVq1a1tSQAAFCM2WXOyKFDh9SiRQsFBQVp6tSpunDhgmVZdHS0ypcvbwkiktSyZUs5OjrqyJEj9igHAAAUYTaPjNxK69at1b59e9WoUUPx8fGaN2+eBg8erI0bN6pEiRI6d+6cKlasaF2Ek5NcXV2VmJhY0OUAAIAirsDDSKdOnSy/Z01gbdeunWW05G5Se0Huy2JHFV4dAAAUZ3Y/tbdmzZq6//77FRsbK0mqXLmyzp8/b7VORkaGLl26lOs8EwAAcPeyexj5+++/dfHiRUvQ8PX1VXJyso4ePWpZ58CBA8rMzJS3t7e9ywEAAEWMzYdpUlJSFBcXZ7l96tQpxcTEyNXVVa6urlq8eLGCgoJUuXJlxcfH66233lLt2rXVunVrSVL9+vXVunVrTZkyRdOmTVN6errCwsLUqVMnzqQBAOAeZHMYOXr0qAYMGGC5nXU9keDgYL3++uv67bfftHXrVl2+fFlVqlRRq1atNGrUKKtrjcyZM0dhYWF67rnn5OjoqCeeeEKTJ08ugO4AAIDixuYw4ufnp+PHj+e6fPXq1be8jwoVKnCBMwAAIInvpgEAAAYjjAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAzlZHQBsIPXg3Np31K4dQAAcBsYGQEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKJvDyOHDhzV06FD5+/vLZDJp165dlmXp6el666231KVLFzVu3Fj+/v569dVXlZCQYHUfgYGBMplMVj8rV668894AAIBix8nWDVJTU2UymdSjRw+FhoZaLfvnn3907NgxDRs2TB4eHkpOTtaMGTM0bNgwbd682WrdkSNHqnfv3pbbLi4u+ewCAAAozmwOIwEBAQoICMhxWbly5bRmzRqrtilTpqhXr146c+aMqlevbml3cXGRm5ubrbsHAAB3GbvPGbly5YocHBxUvnx5q/Z33nlHfn5+6tatm1atWqWMjAx7lwIAAIogm0dGbHHt2jXNmTNHnTp1UtmyZS3t/fv3V6NGjeTq6qro6GjNmzdPiYmJmjBhgj3LAQAARZDdwkh6erpGjRols9msadOmWS0bNGiQ5XcPDw+VLFlSU6dO1ZgxY1SqVCl7lQQAAIogu4SR9PR0vfzyyzpz5ozWrl1rNSqSEx8fH2VkZOjUqVOqV6+ePUq6K9VekHN7bOGWAQDAHSnwMJIVRGJjY7Vu3Trdf//9t9wmJiZGjo6OqlSpUkGXAwAAijibw0hKSori4uIst0+dOqWYmBi5urrKzc1NI0eO1LFjx7RixQpdv35diYmJkiRXV1eVKlVK0dHR+umnn/Too4/KxcVF0dHRCg8PV9euXeXq6lpwPQMAAMWCzWHk6NGjGjBggOV2eHi4JCk4OFihoaH6+uuvJUlPPfWU1Xbr1q2Tn5+fSpUqpaioKC1evFhpaWmqUaOGBg4caDWPBAAA3DtsDiN+fn46fvx4rsvzWiZJDz/8sCIiImzdLQAAuEvx3TQAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFB2+dZe3KHXg/NYtqXw6gAAoBAwMgIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAoJ6MLuGu9HpxL+5bCrQMAgCKOkREAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEPZHEYOHz6soUOHyt/fXyaTSbt27bJabjabtWDBAvn7+8vb21sDBw7UyZMnrda5ePGixowZoyZNmqhZs2aaOHGiUlJS7qgjAACgeLL5OiOpqakymUzq0aOHQkNDsy1/55139P7772vmzJmqUaOGFixYoJCQEEVFRal06dKSpLFjxyoxMVFr1qxRenq6Jk6cqNdee01z58698x7hznGNFABAIbJ5ZCQgIECjR49W+/btsy0zm81at26dhg0bpnbt2snDw0OzZ8/W2bNnLSMof/zxh/bu3avp06fLx8dHzZo10+TJk7V9+3YlJCTceY8AAECxUqBzRk6dOqXExES1bNnS0lauXDn5+PgoOjpakhQdHa3y5cvLy8vLsk7Lli3l6OioI0eOFGQ5AACgGCjQMJKYmChJqlSpklV7pUqVdO7cOUnSuXPnVLFiRavlTk5OcnV1tWwPAADuHZxNAwAADFWgYcTNzU2SlJSUZNWelJSkypUrS5IqV66s8+fPWy3PyMjQpUuXLNsDAIB7R4GGkRo1asjNzU379++3tF25ckU//fSTfH19JUm+vr5KTk7W0aNHLescOHBAmZmZ8vb2LshyAABAMWDzqb0pKSmKi4uz3D516pRiYmLk6uqq6tWra8CAAVq2bJlq165tObW3SpUqateunSSpfv36at26taZMmaJp06YpPT1dYWFh6tSpk6pWrVpwPQMAAMWCzWHk6NGjGjBggOV2eHi4JCk4OFgzZ87U4MGDdfXqVb322mtKTk5W06ZNtWrVKss1RiRpzpw5CgsL03PPPSdHR0c98cQTmjx5cgF0BwAAFDc2hxE/Pz8dP3481+UODg4aNWqURo0ales6FSpU4AJnAABAEmfTAAAAgxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABjKqaDvMDAwUKdPn87W3rdvX02dOlX9+/fXoUOHrJb16dNHb7zxRkGXAgAAioECDyOffPKJrl+/brn9+++/a9CgQerQoYOlrXfv3ho5cqTl9n333VfQZQAAgGKiwMNIxYoVrW6vXLlStWrV0iOPPGJpK1OmjNzc3Ap61wAAoBiy65yRtLQ0ffbZZ+rRo4ccHBws7ZGRkfLz81Pnzp01d+5cXb161Z5lAACAIqzAR0b+bdeuXbp8+bKCg4MtbZ07d1b16tVVpUoVHT9+XHPmzNGJEye0ePFie5YCAACKKLuGkU2bNqlNmzaqWrWqpa1Pnz6W300mk9zc3DRw4EDFxcWpVq1a9iwHAAAUQXY7THP69Gl999136tmzZ57r+fj4SJJiY2PtVQoAACjC7BZGNm/erEqVKumxxx7Lc72YmBhJYkIrAAD3KLscpsnMzNTmzZvVrVs3OTn93y7i4uIUGRmpgIAAVahQQcePH1d4eLiaN28uDw8Pe5QCAACKOLuEke+++05nzpxRjx49rNpLliyp/fv3a926dUpNTdUDDzygJ554Qi+99JI9ygAAAMWAXcKIv7+/jh8/nq39gQce0AcffGCPXQIAgGKK76YBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEM5GV0AjFF7Qe7LYguvDAAAGBkBAADGIowAAABDEUYAAIChmDOCgvV6cB7LthReHQCAYqPAR0YWLVokk8lk9dOhQwfL8mvXrmnatGny8/OTr6+vRowYoXPnzhV0GQAAoJiwy8jIQw89pDVr1lhulyhRwvL7m2++qT179ujtt99WuXLlFBYWptDQUG3YsMEepQAAgCLOLmGkRIkScnNzy9Z++fJlbdq0SXPmzFGLFi0k3QgnHTt21I8//qjGjRvboxwAAFCE2WUCa2xsrPz9/fX4449rzJgxOnPmjCTp6NGjSk9PV8uWLS3r1q9fX9WrV9ePP/5oj1IAAEARV+AjI97e3goPD1fdunWVmJioJUuW6Nlnn1VkZKTOnTunkiVLqnz58lbbVKpUSYmJiQVdCgAAKAYKPIwEBARYfvfw8JCPj4/atm2rHTt2qEyZMgW9OwAAUMzZ/Toj5cuXV506dRQXF6fKlSsrPT1dycnJVuskJSXlOMcEAADc/ex+nZGUlBTFx8fLzc1Nnp6eKlmypPbv36+goCBJ0p9//qkzZ84YNnk1t+9oiR1VuHUAAHCvKvAwMmvWLLVt21bVq1fX2bNntWjRIjk6Oqpz584qV66cevTooZkzZ8rV1VVly5bV9OnT5evry5k0AADcowo8jPz999965ZVXdPHiRVWsWFFNmzZVRESEKlasKEmaOHGiHB0dNXLkSKWlpcnf319Tp04t6DIAAEAxUeBhZP78+XkuL126tKZOnUoAAQAAkviiPAAAYDDCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBRhBAAAGIowAgAADEUYAQAAhiKMAAAAQzkZXQCKp9oLcm6PLdwyAAB3AUZGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYirNpULS8HpxL+5bC2R4AUOgYGQEAAIYq8JGRFStW6IsvvtCff/6pMmXKyNfXV2PHjlW9evUs6/Tv31+HDh2y2q5Pnz564403CrocAABQxBV4GDl06JCeffZZeXl56fr165o3b55CQkK0fft2OTs7W9br3bu3Ro4cabl93333FXQpAACgGCjwMLJ69Wqr2zNnzlSLFi30yy+/qHnz5pb2MmXKyM3NraB3DwAAihm7zxm5fPmyJMnV1dWqPTIyUn5+furcubPmzp2rq1ev2rsUAABQBNn1bJrMzEy9+eabatKkidzd3S3tnTt3VvXq1VWlShUdP35cc+bM0YkTJ7R48WJ7loMihO+2AQBksWsYmTZtmn7//XetX7/eqr1Pnz6W300mk9zc3DRw4EDFxcWpVq1a9iwJAAAUMXY7TPPGG29o9+7dWrt2rapVq5bnuj4+PpKk2Fg+FwMAcK8p8JERs9mssLAwffnll3r//fdVs2bNW24TExMjSUxoBQDgHlTgYWTatGnatm2bli5dKhcXFyUmJkqSypUrpzJlyiguLk6RkZEKCAhQhQoVdPz4cYWHh6t58+by8PAo6HIAAEARV+Bh5KOPPpJ048Jm/xYeHq7u3burZMmS2r9/v9atW6fU1FQ98MADeuKJJ/TSSy8VdCkAipBcJy2PKtw6ABQ9BR5Gjh8/nufyBx54QB988EFB7xYAABRTfFFebvjCNQAACgVflAcAAAxFGAEAAIYijAAAAEMRRgAAgKEIIwAAwFCEEQAAYCjCCAAAMBTXGQEKCtemQWHjOYe7BCMjAADAUIyM4J6T63ekXMjlU6bEJ00AsCNGRgAAgKEYGUGxk9vIhiTFFl4ZxRLfnAugKGJkBAAAGIowAgAADEUYAQAAhiKMAAAAQxFGAACAoQgjAADAUIQRAABgKMIIAAAwFGEEAAAYiiuwAsUI36uTP3letZerzwKGY2QEAAAYipERwEa5jk7c6R2/nsvoBiMbeSsGj5vNI1pFqHYUsGLwfDUCIyMAAMBQjIwA4NMaAEMxMgIAAAzFyAiA25LnGSmFMPfBbnN1UPhyG4mTGI27RzEyAgAADMXICFCI8hxdKLwy8G+3MV/mrjwbhnlCKEIYGQEAAIYijAAAAEMRRgAAgKGYMwLAWMxdKJ6M/LvxnLnrMDICAAAMRRgBAACGIowAAABDGTZn5MMPP9Tq1auVmJgoDw8PTZkyRd7e3kaVAwC4V9zJnBOjrx57l86XMWRkJCoqSuHh4Ro+fLi2bNkiDw8PhYSEKCkpyYhyAACAgQwZGVmzZo169+6tHj16SJKmTZum3bt3a9OmTXrxxReNKAkAgLub0aM6eSj0MJKWlqZffvlFQ4YMsbQ5OjqqZcuWio6OzrZ+ZmamJOnq1at2qaeeS87tqZnVclmQestt73T7XLe90+2pPX/7vtPtqf2O932n2xfn2vPkegfb38m29tr3nW5P7fnb1pb92yDrfTvrfTwvDmaz2VzgFeQhISFBbdq00YYNG+Tr62tpnz17tg4fPqyPP/7Yav2kpCSdPHmyMEsEAAAFpE6dOqpUqVKe6xT5i565urqqTp06Kl26tBwdOfkHAIDiIDMzU9euXZOrq+st1y30MHL//ferRIkS2SarJiUlqXLlytnWd3JyumWiAgAARU/ZsmVva71CH2ooVaqUHn74Ye3fv9/SlpmZqf3791sdtgEAAPcGQw7TDBo0SOPGjZOnp6e8vb21du1aXb16Vd27dzeiHAAAYCBDwkjHjh11/vx5LVy4UImJiWrYsKFWrVqV42EaAABwdyv0s2kAAAD+7a49PeXDDz9UYGCgvLy81KtXLx05ciTP9Xfs2KEOHTrIy8tLXbp00Z49ewqp0vyzpY8RERHq27evmjdvrubNm2vgwIG3fEyKAlv/jlm2b98uk8mkl156yc4V3hlb+5ecnKxp06bJ399fnp6eCgoKKvLPVVv7+N577ykoKEje3t4KCAjQm2++qWvXrhVStbY7fPiwhg4dKn9/f5lMJu3ateuW2xw8eFDBwcHy9PRU+/bttXnz5kKoNH9s7d8XX3yhQYMG6dFHH1WTJk3Up08f7d27t5CqzZ/8/A2z/PDDD2rUqJGeeuopO1Z45/LTx7S0NM2fP19t27aVp6enAgMD9cknn9ilvrsyjNh6ufn//ve/GjNmjHr27KmtW7fq8ccf1/Dhw/Xbb78VcuW3z9Y+Hjx4UJ06ddK6deu0YcMGPfDAA3r++eeVkJBQyJXfvvx+bcCpU6c0a9YsNWvWrJAqzR9b+5eWlqZBgwbp9OnTWrBggXbu3KmwsDBVrVq1kCu/fbb2MTIyUnPnzlVoaKiioqI0Y8YMRUVFad68eYVc+e1LTU2VyWTS1KlTb2v9+Ph4DRkyRH5+fvr000/13HPPafLkyUX2DdvW/h0+fFgtW7bUypUrtXnzZvn5+WnYsGE6duyYnSvNP1v7mCU5OVnjxo1TixYt7FRZwclPH0eNGqX9+/drxowZ2rlzp+bOnau6devap0DzXahnz57madOmWW5fv37d7O/vb16xYkWO648aNcr84osvWrX16tXLPGXKFLvWeSds7ePNMjIyzL6+vuYtW7bYqcI7l58+ZmRkmPv06WOOiIgwjxs3zjxs2LDCKDVfbO3f+vXrzY8//rg5LS2tsEq8Y7b2cdq0aeYBAwZYtYWHh5uffvppu9ZZUNzd3c1ffvllnuvMnj3b3KlTJ6u2l19+2fz888/bs7QCcTv9y0nHjh3NixYtskNFBc+WPr788svm+fPnmxcuXGju2rWrnSsrOLfTxz179pibNm1qvnDhQqHUdNeNjGRdbr5ly5aWtrwuNy9JP/74Y7Zk6+/vrx9//NGepeZbfvp4s6tXryojI+O2LkZjhPz2ccmSJapUqZJ69epVGGXmW3769/XXX6tx48Z644031LJlS3Xu3FnLly/X9evXC6tsm+Snj76+vvrll18sh3Li4+O1Z88eBQQEFErNhaG4vd7cqczMTKWkpKhChQpGl1KgNm3apPj4eIWGhhpdil18/fXX8vT01KpVq9S6dWsFBQVp1qxZ+ueff+yyvyJ/BVZbXbhwQdevX892obRKlSrpzz//zHGbc+fOZTuTp1KlSjp37pzd6rwT+enjzebMmaMqVapYvVEUJfnp4/fff69PPvlEW7duLYQK70x++hcfH68DBw6oS5cuWrlypeLi4jRt2jRlZGQUyRfE/PSxS5cuunDhgvr27Suz2ayMjAw9/fTTGjp0aGGUXChyer2pXLmyrly5on/++UdlypQxqDL7WL16tVJTU/Xkk08aXUqBOXnypObOnasPP/xQTk533duopBuvNz/88INKly6tJUuW6MKFC5o2bZouXryo8PDwAt/fXTcygltbuXKloqKitHjxYpUuXdrocgrElStX9OqrryosLEwVK1Y0uhy7MJvNqlSpksLCwuTp6amOHTtq6NCh2rBhg9GlFZiDBw9qxYoVmjp1qjZv3qzFixdrz549WrJkidGlIR8iIyO1ZMkSvf3223fNlbSvX7+uMWPGaMSIEfabP1EEmM1mOTg4aM6cOZbJ5OPHj9eWLVvsMjpy10U6Wy83L934VHLzKEhe6xstP33Msnr1aq1cuVJr1qyRh4eHPcu8I7b2MT4+XqdPn9awYcMsbVnfFNmoUSPt3LlTtWrVsm/RNsjP39DNzU1OTk4qUaKEpa1evXpKTExUWlqaSpUqZdeabZWfPi5YsEBdu3a1HGYzmUxKTU3Va6+9pmHDht0V30+V0+vNuXPnVLZs2btqVGT79u2aPHmyFixYUGRHYPMjJSVFR48eVUxMjMLCwiTdeK0xm81q1KiRVq9eXSwmtN6Km5ubqlatqnLlylna6tevL7PZrL///lt16tQp0P0V///sm+TncvONGzfWgQMHrNq+++47NW7c2J6l5lt+L6n/zjvvaOnSpVq1apW8vLwKo9R8s7WP9erVU2RkpLZu3Wr5CQwMlJ+fn7Zu3apq1fL46mwD5Odv2KRJE8XFxVl9HffJkyfl5uZW5IKIlL8+/vPPP9kCR1b4Mt8ll0Qqbq83+bFt2zZNmDBBc+fO1WOPPWZ0OQWqbNmy2V5rnn76adWtW1dbt26Vj4+P0SUWiCZNmujs2bNKSUmxtJ04cUKOjo52eT2960ZGpFtfbv7VV19V1apVNWbMGEnSgAED1L9/f7377rsKCAhQVFSUjh49qjfeeMPIbuTJ1j6uXLlSCxcu1Ny5c/Xggw8qMTFRkuTs7CwXFxfD+pEXW/pYunRpubu7W21fvnx5ScrWXlTY+jd85pln9MEHH2jGjBnq16+fYmNjtWLFCvXv39/IbuTJ1j62bdtWa9asUaNGjeTt7a24uDgtWLBAbdu2tRoRKkpSUlIUFxdnuX3q1CnFxMTI1dVV1atX19y5c5WQkKDZs2dLkp5++ml9+OGHmj17tnr06KEDBw5ox44dWrFihVFdyJOt/YuMjNT48eM1ceJE+fj4WF5rypQpY/UpuyixpY+Ojo7ZXlMqVaqU42tQUWLr37Fz585aunSpJkyYoJEjR+rChQt666231KNHD7uM4N2VYeRWl5v/66+/rD59NWnSRHPmzNHbb7+tefPmqU6dOlqyZEmRfmLZ2scNGzYoPT1dI0eOtLqf0NBQjRgxolBrv1229rG4sbV/DzzwgFavXq3w8HB17dpVVatW1YABAzR48GCjunBLtvZx2LBhcnBw0Ntvv62EhARVrFhRbdu21ejRo43qwi0dPXpUAwYMsNzOmtwXHBysmTNnKjExUX/99Zdlec2aNbVixQqFh4dr3bp1qlatmqZPn67WrVsXeu23w9b+RUREKCMjQ2+88YbVB7qs9YsiW/tYHNnaRxcXF7377ruaPn26evTooQoVKujJJ5/Uyy+/bJf6uBw8AAAwVPH9WAkAAO4KhBEAAGAowggAADAUYQQAABiKMAIAAAxFGAEAAIYijAAAAEMRRgAAuEcdPnxYQ4cOlb+/v0wmk3bt2mXzfZjNZq1evVpBQUHy9PRU69attWzZMpvu4668AisAALi11NRUmUwm9ejRQ6Ghofm6jxkzZmjfvn169dVX5e7urkuXLunSpUs23QdhBACAe1RAQIACAgJyXZ6Wlqb58+dr27Ztunz5sh566CGNHTtWfn5+kqQ//vhDH330kSIjI1WvXj1JN77ywFYcpgEAADl64403FB0drfnz5+uzzz5Thw4d9MILL+jkyZOSpK+//lo1atTQ7t27FRgYqMDAQE2aNEkXL160aT+EEQAAkM2ZM2e0efNmLViwQM2aNVOtWrUUEhKipk2bavPmzZKk+Ph4nTlzRjt37tTs2bMVHh6uX375JduXst4Kh2kAAEA2v/32m65fv64OHTpYtaelpalChQqSbkxeTUtL06xZs1S3bl1JN+aQdO/eXX/++afl0M2tEEYAAEA2qampKlGihDZt2qQSJUpYLXN2dpYkubm5ycnJyRJEJKl+/fqSpL/++oswAgAA8q9hw4a6fv26zp8/r2bNmuW4TpMmTZSRkaG4uDjVqlVLkizzSapXr37b+2LOCAAA96iUlBTFxMQoJiZGknTq1CnFxMTozJkzqlu3rrp06aJXX31VX3zxheLj43XkyBGtWLFCu3fvliS1bNlSDz/8sCZOnKhjx47p6NGjeu2119SqVSur0ZJbcTCbzWZ7dBAAABRtBw8e1IABA7K1BwcHa+bMmUpPT9eyZcu0detWnT17VhUqVFDjxo01YsQImUwmSVJCQoKmT5+uffv2ydnZWW3atNG4ceMs80puB2EEAAAYisM0AADAUIQRAABgKMIIAAAwFGEEAAAYijACAAAMRRgBAACGIowAAABDEUYAAIChCCMAAMBQhBEAAGAowggAADDU/wfGMtKUUGy0vgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot histagram of visible_pts_cnt\n",
    "plt.figure()\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "\n",
    "bins = np.linspace(0, 1.6e6, 30)\n",
    "plt.hist([visible_pts_cnt, visible_pts_cnt_single], bins, label=['LoD', 'Single LoD2'], color=['dodgerblue', 'coral'])\n",
    "plt.title(\"Number Distribution of Visible Gaussians\")\n",
    "plt.legend(loc='upper right')\n",
    "plt.grid(False)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gaussian_splatting",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
